


* Title: Behind political affiliation: how moral values, identity politics, and party loyalty have affected COVID-19 vaccination


*************************************************************************************************
***    PRELIMINARY SETTING   ***
*************************************************************************************************
* Clear any data currently in Stata's memory 
 clear all

version 17

* Tell Stata to run the entire do-file without pausing
set more off


*Master
//cap ssc install distplot
//cap ssc install parmest
//cap ssc install cltest


*********************************************************************************
**#  4. Baseline analysis 
*********************************************************************************

**# Table 2.
// Political partisanship and COVID-19 vaccination rates.
* In addition to Table 2, this portion of the code also computes the sharpened q-values for the regression models testing the three hypothesis (i.e., Table 2, Table 4, Table 5)

 
qui{
. use ".\anti-vaxxer dataset - temporary.dta", clear

xtset num_fips date
 global predictors Unemployment_rate_2021 Median_HHIncome2020  poverty_rate prop_lessthan_HSedu pop_density internet_subsc ep_age65 ep_minrty health_insurance_rate staffedallbedsper1000people


* Regression 1. 
  eststo clear
 eststo baseline1: reghdfe series_complete_pop_pct i.rep_majority_2020##i.qdate L7.(total_cases_rate total_deaths_rate) $predictors, a(mdate state i.state_ID#i.mdate) cluster(i.state_ID i.wdate) 
estadd local controlsvar "Yes"
 mat OUTPUT = r(table)
* Display the matrix to check its structure
mat list OUTPUT
* Extract the row named "pvalue" for the i.rep_majority_2020#i.qdate coefficients
matrix pvals_21 = OUTPUT["pvalue", "1.rep_majority_2020#245.qdate".."1.rep_majority_2020#249.qdate"] // alternatively, OUTPUT["pvalue",16..20]
* Display the extracted p-values
matrix list pvals_21



* Regression 2. 
 eststo baseline2: reghdfe series_complete_pop_pct i.partisanship_2020##i.qdate L7.(total_cases_rate total_deaths_rate) $predictors, a(mdate state i.state_ID#i.mdate) cluster(i.state_ID i.wdate) 
estadd local controlsvar "Yes"
mat OUTPUT = r(table)
mat list OUTPUT
matrix pvals_22 = OUTPUT["pvalue", "4.partisanship_2020#245.qdate".."4.partisanship_2020#249.qdate"] 
matrix list pvals_22


 
* Regression 3.
 eststo baseline3: reghdfe series_complete_pop_pct i.rooted_partisanship##i.qdate L7.(total_cases_rate total_deaths_rate) $predictors, a(mdate state i.state_ID#i.mdate) cluster(i.state_ID i.wdate) 
estadd local controlsvar "Yes"
mat OUTPUT = r(table)
mat list OUTPUT
matrix pvals_23 = OUTPUT["pvalue", "3.rooted_partisanship#245.qdate".."3.rooted_partisanship#249.qdate"] 
matrix list pvals_23



* Regression 4.
 eststo baseline4: reghdfe series_complete_pop_pct i.trumpist_county##i.qdate L7.(total_cases_rate total_deaths_rate) $predictors if rep_majority_2020==1, a(mdate state i.state_ID#i.mdate) cluster(i.state_ID i.wdate) 
estadd local controlsvar "Yes"
mat OUTPUT = r(table)
mat list OUTPUT
matrix pvals_24 = OUTPUT["pvalue", "1.trumpist_county#245.qdate".."1.trumpist_county#249.qdate"] 
matrix list pvals_24


* Join the matrices 
matrix pvals_2 = pvals_21 , pvals_22 , pvals_23 , pvals_24
matrix list pvals_2
}
 
 

* Table 2
 esttab baseline1 baseline2 baseline3 baseline4 using ".\Table2.tex", stats(N r2 controlsvar,  labels("Observations" "R-squared" "Health and economic controls")) compress replace b(3) se(3) nonotes nogaps  starlevels(* 0.10 ** 0.05 *** 0.01) mtitle("Model (1)" "Model (2)" "Model (3)" "Model (4)") drop($predictors 2**.qdate ***244.qdate 0.rep_majority_2020*** 1.partisanship_2020*** 1.rooted_partisanship*** 0.trumpist_county***) ///
 title("Political partisanship and Covid-19 vaccination rates") ///
 addnotes("The dependent variable is Fully vaccinated people (\%). Models 1-4 include state, month and state*month fixed effects. Model 4 is conducted only on counties where Republicans obtained the majority of votes in the 2020 presidential election. Health and economic controls include: unemployment rate, poverty rate, median income, proportion of people with less than HS diploma, population density, proportion of people with internet subscription, proportion of people older than 65 years old, proportion of minority, proportion of people with health insurance, health insurance rate number of hospital beds of all types per 1000 people. Q1, Q2, Q3, and Q4 refer to the periods: Jan 1–Mar 31, Apr 1–Jun 30, Jul 1–Sep 30, and Oct 1–Dec 31, respectively. All models use two-way robust standard errors clustered by state and week. Robust standard errors in parentheses. *** p<0.01, ** p<0.05, * p<0.1.") label // drop oncell to have SE below coeff

 

 
***********************************************************************************
**# 5. The role of moral values and partisanship 
***********************************************************************************

**# S2 Table.
. use ".\anti-vaxxer dataset - temporary.dta", replace

* Generating dummies for high and low values
 // drop vhigh* vlow* high*
 foreach var of varlist social_capita_index relative_comm_1518 {

	gen vhigh`var'= 0
	gen vlow`var'= 0
	gen high`var'= 0
	
	qui sum `var' if date1 == "31may2022", detail 
	
	replace vhigh`var'= 1 if `var' >= `r(p75)' 
	replace vlow`var'= 1 if `var' <= `r(p25)'
	replace high`var'= 1 if `var' >= `r(p50)'
} 
 
  label var vhighsocial_capita_index "Social Capital Index, upper quartile"
 label var vhighrelative_comm_1518 "Relative communal values, upper quartile"
 xtset num_fips date
 
 save ".\anti-vaxxer dataset - temporary.dta", replace


  eststo clear
eststo baseline1: reghdfe series_complete_pop_pct c.relative_comm_1518##i.qdate L7.(total_cases_rate total_deaths_rate), a(mdate fips i.state_ID#i.mdate) cluster(i.state_ID i.wdate) 
estadd local county_fe "Yes"
estadd local state_fe "No"
estadd local month_fe "Yes"
estadd local statexmonth_fe "Yes"
estadd local controlsvar "No"


eststo baseline2: reghdfe series_complete_pop_pct i.vhighsocial_capita_index##i.qdate L7.(total_cases_rate total_deaths_rate), a(mdate fips i.state_ID#i.mdate) cluster(i.state_ID i.wdate) 
estadd local county_fe "Yes"
estadd local state_fe "No"
estadd local month_fe "Yes"
estadd local statexmonth_fe "Yes"
estadd local controlsvar "No"


eststo baseline3: reghdfe  series_complete_pop_pct i.vhighsocial_capita_index##i.qdate c.relative_comm_1518##i.qdate L7.(total_cases_rate total_deaths_rate), a(mdate fips i.state_ID#i.mdate) cluster(i.state_ID i.wdate) 
estadd local county_fe "Yes"
estadd local state_fe "No"
estadd local month_fe "Yes"
estadd local statexmonth_fe "Yes"
estadd local controlsvar "No"


 esttab baseline1 baseline2 baseline3 using ".\S2_table.tex", stats(N r2,  labels("Observations" "R-squared")) compress replace b(3) se(3) nogaps   starlevels(* 0.10 ** 0.05 *** 0.01) drop( 0.vhighsocial_capita_index*** 2**.qdate ***244.qdate*** 1.vhighsocial_capita_index relative_comm_1518) ///
 addnotes("The dependent variable is Fully vaccinated people (\%). All three models use with two-way cluster robust standard errors (state and week). Robust standard errors in parentheses. *** p<0.01, ** p<0.05, * p<0.1.") label
 
 

*------------------------------------------------------------------------------
**# Table 3. 
// Impact of political partisanship on vaccination rate for communities with low levels of universalist values and social capital.

. use ".\anti-vaxxer dataset - temporary.dta", replace

 gen lowsocial_capita_index = 1 if highsocial_capita_index == 0
 label var lowsocial_capita_index "Social capital below the mean"
 
  global predictors Unemployment_rate_2021 Median_HHIncome2020  poverty_rate prop_lessthan_HSedu pop_density internet_subsc ep_age65 ep_minrty health_insurance_rate staffedallbedsper1000people

 
  eststo clear
eststo baseline1: reghdfe series_complete_pop_pct i.rooted_partisanship##i.qdate  L7.(total_cases_rate total_deaths_rate) $predictors if vhighrelative_comm_1518 == 1, a(mdate state i.state_ID#i.mdate) cluster(i.state_ID i.wdate) 
estadd local county_fe "No"
estadd local month_fe "Yes"
estadd local state_fe "Yes"
estadd local statexmonth_fe "Yes"
estadd local controlsvar "Yes"
mat OUTPUT = r(table)
mat list OUTPUT
matrix pvals_41 = OUTPUT["pvalue", "3.rooted_partisanship#245.qdate".."3.rooted_partisanship#249.qdate"] 
matrix list pvals_41


 eststo baseline2: reghdfe series_complete_pop_pct i.rooted_partisanship##i.qdate  L7.(total_cases_rate total_deaths_rate) $predictors if lowsocial_capita_index == 1, a(mdate state i.state_ID#i.mdate) cluster(i.state_ID i.wdate) 
estadd local county_fe "No"
estadd local month_fe "Yes"
estadd local state_fe "Yes"
estadd local statexmonth_fe "Yes"
estadd local controlsvar "Yes"
mat OUTPUT = r(table)
mat list OUTPUT
matrix pvals_42 = OUTPUT["pvalue", "3.rooted_partisanship#245.qdate".."3.rooted_partisanship#249.qdate"] 
matrix list pvals_42

 
eststo baseline3: reghdfe series_complete_pop_pct i.trumpist_county##i.qdate  L7.(total_cases_rate total_deaths_rate) $predictors if vhighrelative_comm_1518 == 1 &  rep_majority_2020==1, a(mdate state i.state_ID#i.mdate) cluster(i.state_ID i.wdate) 
estadd local county_fe "No"
estadd local month_fe "Yes"
estadd local state_fe "Yes"
estadd local statexmonth_fe "Yes"
estadd local controlsvar "Yes"
mat OUTPUT = r(table)
mat list OUTPUT
matrix pvals_43 = OUTPUT["pvalue", "1.trumpist_county#245.qdate".."1.trumpist_county#249.qdate"] 
matrix list pvals_43


 eststo baseline4: reghdfe series_complete_pop_pct i.trumpist_county##i.qdate  L7.(total_cases_rate total_deaths_rate) $predictors if lowsocial_capita_index == 1 &  rep_majority_2020==1, a(mdate state i.state_ID#i.mdate) cluster(i.state_ID i.wdate) 
estadd local county_fe "No"
estadd local month_fe "Yes"
estadd local state_fe "Yes"
estadd local statexmonth_fe "Yes"
estadd local controlsvar "Yes"
mat OUTPUT = r(table)
mat list OUTPUT
matrix pvals_44 = OUTPUT["pvalue", "1.trumpist_county#245.qdate".."1.trumpist_county#249.qdate"] 
matrix list pvals_44


* Join the matrices 
matrix pvals_4 = pvals_41 , pvals_42 , pvals_43 , pvals_44
matrix list pvals_4

 esttab baseline1 baseline2 baseline3 baseline4 using "./Table3.tex", stats(N r2, labels("Observations" "R-squared")) compress replace b(3) se(3) nogaps starlevels(* 0.10 ** 0.05 *** 0.01)   ///
 addnotes("The dependent variable is Fully vaccinated people (%). Models 1-4 include state, month and state*month fixed effects. Models 3 and 4 are conducted only on counties where Republicans obtained the majority of votes in the 2020 presidential election. All models include the following health and economic controls: unemployment rate, poverty rate, median income, proportion of people with less than HS diploma, population density, proportion of people with internet subscription, proportion of people older than 65 years old, proportion of minority, proportion of people with health insurance, health insurance rate number of hospital beds of all types per 1000 people. Q1, Q2, Q3, and Q4 refer to the periods: Jan 1–Mar 31, Apr 1–Jun 30, Jul 1–Sep 30, and Oct 1–Dec 31, respectively. All models use two-way robust standard errors clustered by state and week. Robust standard errors in parentheses. *** p<0.01, ** p<0.05, * p<0.1.") label 
  

  
*------------------------------------------------------------------------------  
**# Figure 3. 
// Effect of rooted partisanship and Trumpism on vaccination rate compliance.

. use ".\anti-vaxxer dataset - temporary.dta", replace

xtset num_fips date

*** Full sample
 reghdfe series_complete_pop_pct i.rooted_partisanship##i.qdate L7.(total_cases_rate total_deaths_rate) $predictors, a(mdate state i.state_ID#i.mdate) cluster(i.state_ID i.wdate) 
estimates store full1

 reghdfe series_complete_pop_pct i.trumpist_county##i.qdate L7.(total_cases_rate total_deaths_rate) $predictors if rep_majority_2020==1, a(mdate state i.state_ID#i.mdate) cluster(i.state_ID i.wdate) 
estimates store full2
 
 *** Subsample
  reghdfe series_complete_pop_pct i.rooted_partisanship##i.qdate  L7.(total_cases_rate total_deaths_rate) $predictors if vhighrelative_comm_1518 == 1, a(mdate state i.state_ID#i.mdate) cluster(i.state_ID i.wdate) 
 estimates store sub1

reghdfe series_complete_pop_pct i.trumpist_county##i.qdate  L7.(total_cases_rate total_deaths_rate) $predictors if vhighrelative_comm_1518 == 1 &  rep_majority_2020==1, a(mdate state i.state_ID#i.mdate) cluster(i.state_ID i.wdate) 
 estimates store sub2
 
coefplot (full1, label("")) (full2, label(: Full sample)) (sub1, label("")) (sub2, label(: Counties with high communal value)), ///
keep(*.qdate) drop(2.rooted_partisanship*) xline(0) ///
rename(1.trumpist_county#245.qdate="High-Trump-Support x Q2 2021" 1.trumpist_county#246.qdate="High-Trump-Support x Q3 2021" 1.trumpist_county#247.qdate="High-Trump-Support x Q4 2021" 1.trumpist_county#248.qdate="High-Trump-Support x Q1 2022" 1.trumpist_county#249.qdate="High-Trump-Support x Q2 2022" ///
3.rooted_partisanship#245.qdate="Safe Rep x Q2 2021" 3.rooted_partisanship#246.qdate="Safe Rep x Q3 2021" 3.rooted_partisanship#247.qdate="Safe Rep x Q4 2021" 3.rooted_partisanship#248.qdate="Safe Rep x Q1 2022" 3.rooted_partisanship#249.qdate="Safe Rep x Q2 2022") ///
scheme(s1color) grid(none) ///
xli(-15 -10 -5, lc(gs12) lw(thin)) xtitle("Regression estimated impact in percentage points")


  graph export ".\Figure 3.png", as(png)  replace // name("Graph")
  graph export ".\Figure 3.eps", as(eps)  replace // name("Graph")

 ***********************************************************************************
**# 6. Co-partisanship and vaccination compliance
***********************************************************************************
 use ".\anti-vaxxer dataset - temporary.dta", clear

xtset num_fips date

gen copartisanship_dem =0
replace copartisanship_dem = 1 if rooted_partisanship2000_20==1 & governor_party_num==0

sum relative_comm_1518, detail
gen high_comm1518=0
replace high_comm1518=1 if relative_comm_1518 > 0.6

 global predictors Unemployment_rate_2021 Median_HHIncome2020  poverty_rate prop_lessthan_HSedu pop_density internet_subsc ep_age65 ep_minrty health_insurance_rate staffedallbedsper1000people
 
 
 
 
*-------------------------------------------------------------------------------
**# Table 4.
// Co-partisanship, communal values, and vaccination compliance.
 
* Regression 1. 
tab co_partisanship rooted_partisanship2000_20

  eststo clear
 eststo regression1: reghdfe series_complete_pop_pct i.co_partisanship L7.(total_cases_rate total_deaths_rate) $predictors if rooted_partisanship2000_20==3 , a(mdate i.state_ID) cluster(i.state_ID i.wdate) 
 estadd local controlsvar "Yes"
estadd local subsample "Rooted Rep"
mat OUTPUT = r(table)
mat list OUTPUT
matrix pvals_51 = OUTPUT["pvalue", "1.co_partisanship"] 
matrix list pvals_51


tab copartisanship_dem rooted_partisanship2000_20

* Regression 2. 
 eststo regression2: reghdfe series_complete_pop_pct i.co_partisanship L7.(total_cases_rate total_deaths_rate) $predictors if rooted_partisanship2000_20==1, a(mdate  i.state_ID) cluster(i.state_ID i.wdate) 
 estadd local controlsvar "Yes"
estadd local subsample "Rooted Dem"
mat OUTPUT = r(table)
mat list OUTPUT
matrix pvals_52 = OUTPUT["pvalue", "1.co_partisanship"] 
matrix list pvals_52



* -------------------------------------------------------------------------------
***** Co-partisanship, compliance and moral values 

* Regression 3. 
 tab co_partisanship rooted_partisanship2000_20

 eststo regression3: reghdfe series_complete_pop_pct i.co_partisanship##c.values_communal1518_shrink L7.(total_cases_rate total_deaths_rate) $predictors if rooted_partisanship2000_20==3 , a(mdate i.state_ID) cluster(i.state_ID i.wdate) 
 estadd local controlsvar "Yes"
estadd local subsample "Rooted Rep"
mat OUTPUT = r(table)
mat list OUTPUT
matrix pvals_53 = OUTPUT["pvalue", "1.co_partisanship"], OUTPUT[ "pvalue","1.co_partisanship#c.values_communal1518_shrink"]  
matrix list pvals_53

* Regression 4. 
 eststo regression4: reghdfe series_complete_pop_pct i.co_partisanship##c.values_communal1518_shrink L7.(total_cases_rate total_deaths_rate) $predictors if rooted_partisanship2000_20==1, a(mdate i.state_ID) cluster(i.state_ID i.wdate) 
 estadd local controlsvar "Yes"
estadd local subsample "Rooted Dem"
mat OUTPUT = r(table)
mat list OUTPUT
matrix pvals_54 = OUTPUT["pvalue", "1.co_partisanship"], OUTPUT[ "pvalue","1.co_partisanship#c.values_communal1518_shrink"]  
matrix list pvals_54


* Join the matrices 
matrix pvals_5 = pvals_51 , pvals_52 , pvals_53 , pvals_54
matrix list pvals_5


 esttab regression1 regression3 regression2 regression4 using ".\Table4.tex", stats(N r2 controlsvar subsample, labels("Observations" "R-squared" "Health and economic controls" "Subsample")) compress replace b(3) se(3) nonotes nogaps  starlevels(* 0.10 ** 0.05 *** 0.01) mtitle("Model (1)" "Model (2)" "Model (3)" "Model (4)")  drop(0.co_partisanship*** $predictors) ///
 addnotes("The dependent variable is Fully vaccinated people (%). Models 1-4 include state and month fixed effects. All models use two-way robust standard errors clustered by state and week. Models 1 and 2 are runned over a sub-sample of 2000-20 stronghold Republican counties, whereas models 3 and 4 are runned over a sub-sample of stronghold Democrat counties. Communal values referes to Abs. importance of communal moral values (2015-2018). Health and economic controls include: unemployment rate, poverty rate, median income, proportion of people with less than HS diploma, population density, proportion of people with internet subscription, proportion of people older than 65 years old, proportion of minority, proportion of people with health insurance, health insurance rate number of hospital beds of all types per 1000 people. Robust standard errors in parentheses. *** p<0.01, ** p<0.05, * p<0.1.") label






 
* -------------------------------------------------------------------------------
**# Figure 4. 
// Co-partisanship, vaccination compliance and moral values 

 use ".\anti-vaxxer dataset - temporary.dta", clear

 xtset num_fips date
 
 * Regression 1. 
 qui{ 
 	eststo clear
 eststo regression2: reghdfe series_complete_pop_pct i.co_partisanship L7.(total_cases_rate total_deaths_rate) $predictors if rooted_partisanship2000_20==1, a(mdate  i.state_ID) cluster(i.state_ID i.wdate) 

  margins co_partisanship, saving(fil1, replace)
  
* Regression 2. 
 eststo regression1: reghdfe series_complete_pop_pct i.co_partisanship L7.(total_cases_rate total_deaths_rate) $predictors if rooted_partisanship2000_20==3 , a(mdate i.state_ID) cluster(i.state_ID ) // cluster only by i.state_ID to get SE when computing margins

 margins co_partisanship, saving(fil2, replace) 

 
combomarginsplot fil1 fil2 , labels("Dem counties" "Rep counties") scheme(s1color) plot1opts(lcolor(dknavy) mcolor(dknavy)) ci1opts(color(dknavy))  plot2opts(lcolor(cranberry) mcolor(cranberry))  ci2opts(color(cranberry)) xtitle("=1 if co-partisan state governor") title("Panel a: predictive margins of co-partisanship with 95% CIs") yli(30 35 40 45 50,  lc(gs12) lw(thin)) ytitle("Effects on vaccination rate") // for the graph see: https://stats.oarc.ucla.edu/stata/faq/how-can-i-understand-a-categorical-by-continuous-interaction-stata-12/
 }
 
  graph export ".\Figure 4A.png", as(png)  replace // name("Graph")
  graph export ".\Figure 4A.eps", as(eps)  replace // name("Graph")


* Regression 3. 
qui{
 eststo regression4: reghdfe series_complete_pop_pct i.co_partisanship##c.values_communal1518_shrink L7.(total_cases_rate total_deaths_rate) $predictors if rooted_partisanship2000_20==1, a(mdate i.state_ID) cluster(i.state_ID i.wdate) 

 margins , dydx(co_partisanship) at(values_communal1518_shrink=(-2(.5)2.5)) vsquish  saving(fil3, replace) // noestimcheck

* Regression 4. 
 eststo regression3: reghdfe series_complete_pop_pct i.co_partisanship##c.values_communal1518_shrink L7.(total_cases_rate total_deaths_rate) $predictors if rooted_partisanship2000_20==3 , a(mdate i.state_ID) cluster(i.state_ID i.wdate) 

 margins , dydx(co_partisanship) at(values_communal1518_shrink=(-2(.5)2.5)) vsquish  saving(fil4, replace) // noestimcheck

// marginsplot, recast(line) recastci(rarea) yline(0)


combomarginsplot fil3 fil4 , labels("Dem counties" "Rep counties") scheme(s1color) recast(line) recastci(rarea)  plot1opts(lcolor(dknavy%20) mcolor(dknavy%20)) ci1opts(color(dknavy%20)) plot2opts(lcolor(red%20) mcolor(red%20))  ci2opts(color(red%20)) yline(0) title("Panel b: average marginal effects of co-partisanship with 95% CIs")  ytitle("Effects on vaccination rate")  // recast(line) recastci(rarea)
}
  graph export ".\Figure 4B.png", as(png)  replace 
  graph export ".\Figure 4B.eps", as(eps)  replace 

 
 
 
 
 * -----------------------------------------------------------------------------
**# 5° RC: Testing for FDR
 ***   Export the matrix with the p-values for later use 
 matrix pvals_all = pvals_2 , pvals_4, pvals_5
matrix list pvals_all
matrix pvals_all_t = pvals_all'
 
 putexcel set ".\pvals_all.xlsx", sheet("pvals_all_t") replace
putexcel A1=matrix(pvals_all_t)
 